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^ , 1 Introduction 

>0 , In recent years radial basis functions (RBFs) have become a popular tool for 

scattered data approximation and interpolation and for solving partial differential 
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equations. Two excellent recent books covering practical and theoretical issues are 
Fasshauer and Wendland [18]. A function # : R'^ — >• M is said to be radial if 
there exists a function (p ■ [0,oo) — >■ K such that <?'(x) = (?!)(||x||2) for all x G K''. 
Then we can define an RBF for a given centre x^ G R as 

*i:(x) = (?!)(||x-x,||2). 

For example, an interpolant I for the scattered data interpolation problem, 
where we are given data (pi.j,yj),j = l,...,n, with Xj e R and j/j € R, can be 
constructed as the linear combination 

n 

2(x) = ^c,<?,(x), xeR^ (1) 

where the coefficients ci , . . . , Cn are chosen so that 

A^j)=Vj^ j = ^,---,ri. (2) 

If <? is positive definite, then ([1]) and ([2]) have a unique solution. We recall 
that a continuous function / : R — >■ R is positive definite (some would say strictly 
positive definite) if for any n distinct points xi, . . . , Xn £ R , the quadratic form 

n n 

^^a,aj/(x, -Xj) 
i=ij=i 

is positive for all a = [ai, . . . , a„]^ G R"\{0}. 

Most geophysical applications are in R'^, but in other areas, such as machine 
learning, different spatial dimensions may occur (see, for example, [3]). 

Broadly speaking RBFs can be categorised as either globally supported or com- 
pactly supported. The first category includes Gaussians and multiquadrics [6] . Both 
have a scale parameter (also known as a shape or tension parameter), the selection 
of which is still a major ongoing research topic. 

The second category includes Wendland functions [TT], Buhmann RBFs [1], 
Wu's RBFs [H] and "Euclid's hat" [81IT3]. The most widely used in this category 
are the Wendland functions, which are piecewise polynomial compactly supported 
functions. They have the minimum polynomial degree for any level of smooth- 
ness, and have a strictly positive Fourier transform. The Wendland functions were 
originally derived for integer-order Sobolev spaces in odd dimensions in Wendland 
[17| and were then extended to even dimensions in Schaback [ll] (Schaback called 
these the "missing" Wendland functions). They are uniquely defined for a given 
spatial dimension d and a smoothness parameter k (up to a constant multiplier). 
All the Wendland functions are equal to zero outside [0,1]. 

Our purpose in this paper is to consider the limiting case of both the orig- 
inal and the missing Wendland functions as the smoothness parameter k goes 
to infinity. In Figure [U we can see the original Wendland functions in R^ for 
k = 1, . . . , 5, where we have normalised the functions to have value 1 at the ori- 
gin. One can clearly see faster decay as k increases, which suggests the need for 



Wendland functions with increasing smoothness converge to a Gaussian 3 

a change of variable when considering the hmit as k approaches infinity. In Sec- 
tion [2] we give some necessary background on radial basis functions and Fourier 
transforms. Then in Section [3] we define the Wendland functions cf)^ fe, where k is 
a non-negative integer for the original Wendland functions and a positive half- 
integer for the missing Wendland functions. We also define normalised Wendland 
functions 4>i,k s-nd normalised equal area Wendland functions tpi,k, where the latter 
are obtained by a linear change of scale in the argument so as to ensure equal 
areas under their graphs for all values of k. In Section H] we consider the limit of 
the normalised equal area Wendland functions as the smoothness parameter goes 
to infinity. Our main theorem is Theorem 14.31 which states that the normalised 
equal area Wendland functions converge uniformly to a Gaussian on the real half- 
line. Section [5] provides numerical evidence to illustrate the theoretical results, as 
well as an example of interpolation using the normalised Wendland functions and 
the normalised equal area Wendland functions. This example compares results 
obtained with the Wendland functions of different smoothness as well as the lim- 
iting Gaussian. The results of the experiment show only small differences in errors 
when using normalised equal area Wendland functions with different smoothness 
parameters, but the condition number of the linear system grows quickly as the 
smoothness parameter increases. A cautious conclusion might be that there is little 
incentive to consider using high values for the smoothness parameter. 

We note that the similarity of the normalised Wendland functions to a Gaussian 
has been mentioned in [11] and of the normalised equal area Wendland functions 
to a Gaussian in [7], in both cases for 'E? with k = 1. No theoretical explanations 
were given for these observations. 




03_l(^) <^,^2^r) 03 3(r) --- 0^^^(r) 0^ ^(r) 



Fig. 1: The (original) Wendland functions (t>3.i{r), 4>i^2{r), 4>B.'i{r), 06, 4('") and 
07,5 (r) normalised to have value 1 at r = 0. 
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2 Background 



In this section, we will provide background material on several types of radial basis 
functions, the Fourier transform and the inverse Fourier transform. 



2.1 Wendland functions 



Wendland functions were originally introduced in [17] and then more cases were 
added in (Tlj. We will refer to the Wendland functions from [T7] as the original 
Wendland functions and the Wendland functions from [14] as the missing Wendland 
functions. A thorough investigation of both types of Wendland functions in terms of 
hypergeometric functions and other special functions is the focus of Hubbert [10] . 
We firstly define the original and missing Wendland functions and then discuss 
the properties of each. 

Definition 2.1 With spatial dimension d, smoothness parameter k > 0, and i : = 
Lf + fcj+l, let 

I , \ , f sil-sYis^^r^f-'^dsfoT 0<r<l, 
for r- > 1. 

This defines the original Wendland functions when k is a positive integer and the 
missing Wendland functions when k is a positive half-integer. 

Both the original Wendland functions and the missing Wendland functions are 
continuous on [0,cxj). Note that this choice of ^ is the minimum value that ensures 
that the resulting functions are positive definite. Since i? = [| + fcj + 1 it follows 
immediately that for fixed d 

£ r^ k as k -^ oo, (4) 

where x ^ y denotes asymptotic equality, that is ^ — )> 1. If we wish to use a differ- 
ent support, this can be easily achieved through scaling the function argument: if 
(/)(r) has support [0, 1] then with the change of variable y = Sr, ipiv) '■= 4> (f ) has 
support [0,5]. 

The original Wendland functions are constructed with the minimal degree poly- 
nomial for a given smoothness that gives rise to a d-dimensional positive definite 
function. They are unique up to a multiplicative constant when fc > 0. Another 
important property is that for an odd space dimension d and for k a non-negative 
integer the function 

7^(x,y) = 0,^fc(||x-y||), x,yGR'', (5) 

is the reproducing kernel of a Hilbert space which is norm equivalent to the Sobolev 



d+l 



space H^-+^{R'') [H Chapter 10]. 
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We give explicit formulae for the original Wendland functions for d = 3 and 
k = 1, 2, . . . , 5 in Table[T] where = denotes equality up to a positive constant factor. 
The support of all the original Wendland functions is [0, 1]. 



k 


Original Wendland function 


1 


<t>3.i{r) = {i - r)X{4r + 1) 


2 


</'4,2(r-) = (l-r)'5^(35r2 + 18r + 3) 


3 


(As.sCr-) = (1 - r)8.(32r3 + 25r-2 + 8r + 1) 


4 


</>6,4(r) = (1 - r)iO(429r4 + 450r^ + 210r2 + 50r + 5) 


5 


^7^5(r) = (1 - r)i2(2048r5 + 2697r4 + 1644r3 + 566r2 + 108r + 9) 



Table 1: The original Wendland functions cj>e,k{'^') for d = 3 and k = 1, 



Schaback [Tl] extended Wendland's original approach to cover the missing 
Wendland functions, which are the reproducing kernels of integer order Sobolev 
spaces in even dimensions d. An important distinction between the original Wend- 
land functions and the missing Wendland functions is that the missing Wend- 
land functions, whilst still being compactly supported, now have logarithmic and 
square-root multipliers of polynomial components. We give explicit formulae for 
the missing Wendland functions for d = 2 and fc = i , | and | in Table [2] The 
support of all the missing Wendland functions is r G [0, 1], while 



L{r) := log 



1 + Vl-r^ 



and S(r) := a/T^ 



It is proved in [T3] that with £ = [^ + k\ +1 and k being any positive half-integer. 



Missing Wendland function 



</,2 1 (r) = 3r^L{r) + (2r^ + l)S{r) 

<t>3,lir) = 15r*(6 + r2)L(r) + {81r* + 28r2 - 4)5(r) 

'I'a 5 (r) = (945r» + 2520r^)L(r) + (256r* + 263r« + egOr* - 136r2 + 16)5(r) 

' 9 



Table 2: The missing Wendland functions 4'i,k{'i') for d = 2 and A; = ^, |, |. 



for an even dimension d the kernel in ([5]) is the reproducing kernel of a Hilbert 
space that is norm equivalent to the Sobolev space 77^""*" (M ). 

Hereafter when we refer to Wendland functions we will mean both the original 
and missing Wendland functions. They will both be denoted by 4>e^k{r) where 



d 



+ k 



+ 1 
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and r G [0, 1], with k an integer for the original Wendland functions and a half- 
integer for the missing Wendland functions. 



2.2 Fourier transforms 

The proofs of the main results rely heavily on Fourier transforms. This section 
provides definitions and outlines some key properties. For further information, we 
refer the interested reader to [151118] . 

With the Fourier transform of / G Li[Mr) defined as 



/(w) := {2^)-i [ /(x) 



e dx, w £ K , 



it is well known that the Fourier transform of a radial function <P G Li (R ) nC(M ) 
is also radial and is given by ^(w) = J"£;(/f)(||w||2) where 



Td(p{z) := z 2 / (l){y)y^ Jd_^{zy)dy, (6) 

Jo ^ 

and Jv{y) denotes the Bessel function of the first kind with order v. 

From the Fourier inversion theorem applied to radial functions, we know that 
if <?e LiiW^) with<5(x) = <;/)( ||x|| 2), (^: [0,oo)^R, and if <? G Li(R'^), then 

d f°° d 

<t>{y) = y^'~ Td4>{z)z^ Jd_^{yz)dz. (7) 



We also recall that if / G Li(R'') n C(R'') is positive definite then its Fourier 
transform is in Li(R ) and is non- negative [181 Corollary 6.12]. 

The Gaussian radial basis function with scale parameter a > 0, which we 
denote by Ga{y), is given by 

and its Fourier transform is given by 

Ga{z) = ^-^e-^, zgR. 
(2a)f 

We will first define the generalised hypergeometric function and then state 
the Fourier transform of the Wendland functions. Further details on generalised 
hypergeometric functions can be found in [1] and [2] . 

Definition 2.2 The generalised hypergeometric function pFq{ai, . . . , ap; hi, ... , bq\ x) 

is 

„/ , 1 \ sr^ {ai)n- ■ ■ {ap)n x^ 
pFq{ai,...,ap;bi,...,bq;x) := > -r— ; r— ; j-, 

^ (bi)n--- (bq)n n\ 
n=Q 

where none of bi, . . . ,bq is a negative integer or zero and where 

{c)„:=c{c+l)---{c + n-l) = £^j±^, n>l (8) 
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denotes the Pochhammer symbol, with (c)o = 1. When p < q the series converges for all 
finite X and defines an entire function. When p = q+ 1 the series converges absolutely 
for \x\ < 1, and also at x = 1 if 

q p 

^ 6j - ^ Oj > 0. 

Theorem 2.3 The d-dimensional Fourier transform of <j)£j. is given by 
-r ^ (^ nii^k r^ (d+l ^, l±d+l l + d + 2^, z^ 

where 



2''+i r{e + i)r (!i±^ + k) 






Proof This result can be found in [20| . however for completeness, we provide it 
here in more detail. With ([3} and (|6]) and with the change of variables s = y/t and 
a; = f in the second line, we have that 



1 1 

, fe-i 



J'd^iM^) = 2k~ir(a) / /(i-*)^ *(*'-?/") y'^|_i(^y)dtdy 



2fc-ir(a) 

y 

11 



fc-1 



b b 
From P 6.567.1], we have that 

[ 3;"+' (l - x^Y J^bx) dx = 2''r(/i + l)6-(''+i^ J^+^+i(6), b > 0, 

which means we can simplify our expression to 

1 
^d<l>i,k{z) = z^~^ \ x^^^(\-xf3^^^^_^{zx)dx. 
b 

With the following identity |9l 6.569] 



(1 — a;)^ Ji/(aa;)da;: 



r(i. + i)r(A + /i + 1/ + 1) 



X2F3(^ 2 ' 2 ''^+^' 2 ' 2 '-T 

the result follows after applying the duplication formula for the gamma function 
[11 6.1.18] 

r{z)r[z-r^^=2^-^^^r{2z), (9) 

and noting that in our case, the ^F-i hypergeometric function simplifies down to a 
\Fi since the first parameters in the numerator and denominator are equal. 
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3 Normalised and normalised equal area Wendland functions 



Hubbert [10] starts with the integral representation IQ to express the Wendland 
functions in terms of Legendre functions. Equation (3.4) in [10] states that for 
re (0,1] 

(10) 
Now we apply the following identity [1] 15.3.4] 

2Fi{a,b;c;z) = (1 - 2)"°2-Fi (a,c~b;c- _ j (11) 

to (|10|1 . which gives us, for r G [0, 1], (since we recover the case of r = by right 
continuity) 

0,,(r)= ,,5/+^^ n-ry^\FJi/-±l;i + k+l;l-A. (12) 

We will need the normalised Wendland functions, which we denote by 4>e^k{i~), 
where 

'^^^'^^^^ ^= 0;:;(o)' ^''^ 

so that <l>i k{0) = 1- 

Theorem 3.1 

- ,, (^§ff^^il-rr+'^.F,U,i±^;e + k+l-l-rA0<r<l, 

[O r > 1. 

(14) 

Proof To calculate 0^ fc(O) we need the value of the hypergeometric function in 
(|12|) at the argument 1 (since it has argument 1 — r^). From [Tj 15.1.20] we have 
the identity 

2Fi(a,b;c;l)= ^i")^j^" .^~"j , c - 6 - a > 0. (15) 

1 (c — 0)1 (c — a) 



Applying ([TS]) to ^^ shows that 
^£,fe(0) 



r(^ + i)r(fc+i) 



2^+^T(| + fc + i)r(| + fc + i)- 

Using (O twice - firstly for r(| +fc+ i)r(| + A:+ 1) and then for r{k)r (fc + i) - 
and with several terms cancelling out in the numerator and denominator, we get 

^ rm r(^ + i)r(2fc) 

'^''''^^> - 2>^~^r{k)r{£ + 2k + iy ^^^> 



Wendland functions with increasing smoothness converge to a Gaussian 
Now using (|12p and (|16p we can see that for r G [0,1] 

which coincides with (|14|) . n 

We will also need the following result for the area under tlie Wendland func- 
tions. 

Lemma 3.2 The area under the graph of (t>£j^(r),r > 0, is 

_ 2'^r{£+i)r{k + i) 
"''' - r{i + 2k + 2) ■ ("^ 

Proof Since the Wendland function is equal to zero for r > 1 we need consider only 
the integral over [0, 1]. Then 

«Afc = / '?^f,/c('")d'- (18) 







n^+1) f n 2.£+fe ,^ fe i + l.„ , ,, , ,., „2 



(1-r^r+SFi 1 ^-^;^ + fc+l;l-r^ dr 



2i+kr{i + k + i) Jq ^ ' -" '\r 2 

■/ J/ (i-y) " 2i^i I 2'^"'^ + '^ + ^'^ )^^ 

h,k, (19) 



2e+k+ir{£ + k + l)Jo 

r{£ + i) 



2^+k+^r{£ + k + l) 



by using (|12p in the second line, the substitution y = 1 — r" in the third line and 
defining Igj^ as the integral in the third line. Using the identity [T] 15.3.3] 

2F1 (a, b; c; 2) = (1 - 2)'"°"'' 2^1 (c - a, c - fe; c; z), (20) 

we can rewrite /^ j, as 

^e,k = l y'+''{l-y)\Fif^ + k + l,^ + k + ^;£ + k + l;y\dy. (21) 

Now we can use the identity 7.512(4)] 

[\--Hi-^r~'.F,ia,p-,r,^)d. = £W£}Ilr+^^ (22) 

Jo -r(7 + p-a)r(7 + p-/3) 

which is valid for Re 7 > 0,Re p > and Re(7 + p-a-^) > 0. With 7 = ^ + ^ + 1, 
p = k + l, a = ^ + k + l and /3 = ^-^ + k and noting that F (^) = i/tt, we get 

, _ r(^ + fc + i)r(fc + i)0F 

'^^'^-r(fc + | + i)r(fc + | + i)- (''^ 

Then using (|9l) and substituting 1^ into ([19]) we get dTT]). D 
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Now we can define what we will call normalised equal area Wendland functions. 
These are normalised Wendland functions with a linear change of variable (which 
depends on k and a positive real constant a) such that for a given a all the 
normalised equal area Wendland functions have area equal to the area under 
exp(— ay^) over the real half-line. The value of a can be chosen for the conve- 
nience of the user. We will denote the normalised equal area Wendland functions 
by V'f.fc- 

Theorem 3.3 With a any positive real constant and 

^^■^^") •= 2^r(fc + i) ' ('^) 

the normalised equal area Wendland functions 



I for y > 5pk{a) 



(25) 



have area under the graph equal to -^JthJoi, the area under exp(— ay ) over the real 
half-line. 



Proof With (|17|) and (|16p . we can see this as follows 

^„ , I 

5l,k{oi) [^ 



't>i,k ( ^ — 77X ) '^V = ^e,k{(^) / </><?, fe(r)dr 



0f,fc('^)dr 



<t>e,kW Jo 
<5^,fc(") 






o-£,k 



2^ 

f°° 2 

/ exp{-ay )dy. 
Jo 



In Figure [2] we plot the normalised equal area Wendland functions tp£ i^ for k 



We emphasise that the normalised equal area Wendland functions satisfy both 
''Pi fe(0) = 1 as well as 

llV'<?,fclli = / e "^ dy, 

■10 

where a can be any constant. 

We will also need the following results. 



Lemma 3.4 For k > min (f , l), 



J,,fe(a) < ^. (26) 

va 
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V3, ,(>■) V4 2(3') V5,3(>'' 1'6,4(J') V7 sl;-) — -e->^ 



Fig. 2: The normalised equal area Wendland functions V'3,i(j/)i V'4,2(y), i^5,3{y), 
V'6,4(y),^7,5(y) with a = 1 and exp(-y^). 



Proof From [16| we have the double inequality 

< 



i-« r{x + s) ^^ 



for < s < 1 and a; > 0. With s = ^ and using r(A: + 1) = kr{k), this gives 

' + 2fc + i)r(fc + i) 



•^^^.fcCQ!) 



2^ r(A:+l) 



(£ + 2fc + 1) 


2Vka 

{[^ + k\ +2k + 2) 



2Vka 



3^/fc 



Lemma 3.5 Let 77 > 0. The function frj : (0, 00) —> R defined by 



Myy-=^' .>o. 



is non-deereasing on (0,cxj). 



(27) 



Proof We provide a proof for the convenience of the reader. Consider Fn{y) 
log/r,(2/). Since f^{y) > and 

dFnJy) _ 1 d/„(y) 
dy /.)(?/) dy 
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we can consider the derivative of Fjj{y) since it will have the same sign as the 
derivative of /,, (y) . Now 

where ^o is the digamma function with the series representation [I] 6.3.16] 

oo 

^o(y + l) = -7+V , ^^ . , 2/ ^-1,-2,-3,... 

n=l ^ "' 

where 7 is the Euler-Mascheroni constant. This gives 



E 



which converges and is positive on (0,cxj), since ry > 0. 



4 Limit of the Wendland functions as fc — >^ cx) 



In this section we derive the limit of the normalised equal area Wendland functions 
as fc — ^ 00 for the original and missing Wendland functions. We start with a 
convergence result for the Fourier transforms. 



Theorem 4.1 Let a be a positive real constant, k be a non-negative integer or positive 
half-integer and ipi,k be the normalised equal area Wendland functions defined by (|25|) 
with (, given by (|2.1|) . Then 

liu, rd^t^k{z) = Gc.{z) (28) 

uniformly for z in an arbitrary bounded subinterval of the positive numbers. 

Proof Firstly we express the Fourier transform of -i/)£ j. in terms of the Fourier 
transform of <^£ j,. Writing 5g j, for Sg ki'^) ^md using the transformation y = rSf i; 
gives, from ^, 

Jo 

= z^~~- Sl.k I 4>e.k{r) {rS(,k)^ J-i_-^{zr5e^h)'ir 
Jo ^ 

= ^^~" ^«,fe' 7; — rrn: {^d4>e,k) {Sg^kz) 
[Se,kz)' 2 

= Si,k {^d4>e,k) {S£,kz)- (29) 
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Next we use Theorem 12.31 and (|16p to write the Fourier transform of ipi_}^ in terms 
of a hypergeometric function as 

5lk2-ir{t + 2k + i)r{k)r{d+2k) 
'^ ' r{2k)r{^ + k)r{£ + d + 2k + i) 

^ , d+l , , £ + d + 2k + l £ + d + 2k + 2 Sl^z-" 
F2\—;T-+k; , ; -7—], (30) 



iJ 



with 5ik = 5£k{a). Then the equivalent series representation is as follows 

_^ , (,x_o-#y- r{d+2k + 2n)r{£ + 2k+l)rik) ^d+2n {~^) 

^ '''^'''^ ^^r{2k)r{£+2k + l + d + 2n)r{k + ^+n) ^^'^ n! 

00 / 2 \ ^'' 

= 2-t^«;„{fc) -^ , (31) 

where 



r{d + 2fc + 27%)r{£ + 2fc + i)r(fc) ^rf+2„ 

r{2k)r{£ + 2fc + 1 + d + 2n)r (fc + 1 + i 



un{k) := ^,„,; :.;,.,,: , . ^ ^ ;, , 7 , ^ ^tX- (32) 



To interchange the limit as fc — >■ cxa and the infinite sum, we need to prove that 
this sum is dominated by an absolutely convergent series (Lebesgue's dominated 
convergence theorem) . Using (|9]) twice together with Lemma 13.41 and the bound 
[51 5.6.8] 

f|^<^, x>0,6-a>l,a>0, (33) 

we have for fc > 1 

^ ^ 2'^+^"r(| + fc + n + l)r(£ + 2fc + l)<5,^l^" 

'^''^ ' r(fc+i)r(£+2fc + l + d + 2n) n\ 

2''+2"r(| + f +n)r(3fc) <^'" 



r(^)r(3fc + d + 2n) 



d+2nj^ ( d I 3A: i ^'l o3/c-l p /'3fc\ p /sfc+l 



2''+^"r § + 



2^2 



+ n)23'=-ir(f)r(^) 5^rf+2 



r (M) 23fc+d+2„-l;. (^3k±d ^ ^^ ^ ^3fc^+l ^ , 



1 /3x^y+'" 



^ 3fc 
^ 2 



gN 2+" X 



Un, 
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where for the inequahties above we used 

r{£ + 2k+l) r{3k) 



r{£ + 2k+l + d + 2n) " r{3k + d + 2n) 
and for A: > 1 

r{l + k + n + ^) ^ r(| + f + n) 
r(fc + i) - r(f) ' 

both of which follow from Lemma [3^ 



/ 2 \ " 

The ratio test shows that ^^^ Un (—^] is absolutely convergent. Therefore 
we can take the liinit as fc — >■ oo inside the infinite sum, giving 



lim J-dipe,kiz) = 2 ^ y lim Wn{k) 
Using the following asymptotic result from [5] 



z^ 



we can see that 



and hence 



r{x + g) a-b 



lim wn{k) = 

fc-i-oo ■nln.2+" 



(34) 



lim Tdi'eA^) = (2a) = V ^ ^ 



n! 

n=0 

= (2a)^2 g-fe 

= Ga(^), 

which proves pointwise convergence. Uniform convergence follows since the inter- 
val is bounded. □ 



Before we consider convergence of tpi,k we establish the following technical lemma. 
Lemma 4.2 For arbitrary a > 0, 

^°° Td^i^z) z"-^ dz = j^ Ga{z) z"-^ dz = 2i-^r [^ , 
where (. is given by (|2.ip . 
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Proof From [121 8.4.48.1], we have the foUowing identity 

r{s)ria^s)r{bi)r{b2) 



pOO 

/ t''"\F2{a; 61,62; - 
Jo 



-t) dt : 



r(a)r(6i-s)r(62-s) 



< Re(s) < min ( -(61 + 62 - a) + -, a ) . (35) 

In our case, (|35[) is vahd at least for k > d since then 

mm (Ul + ^ + k + Ij ,^-^ + k \ >d. 

We can calculate the integral of Tdipe.k over the real half-line as follows. We use 

(|9]) again, and with the change of variable w = -^ — in the second step we see 
that 

' , , d-i^ 5f,2-ir{e + 2k + i)r{k)r{d + 2k) 



r(2fc)r (f + k)r{i + d + 2k + 1) 



P ^ d+l , , l + d + 2k + l l + d + 2k + 2 5lf^z^\ ^_ 

^^ ^^\-^+^-^ — 2 — ' — 2 — ;-^j^ "^^ 

Sjf. 2-ir{£ + 2k + l)r{k)r{d+2k) 2''-i 

r{2k)r{§ + k)r{£ + d + 2k + i) &fj^ 

r T, fd+'^ , , e + d + 2k + l £ + d + 2k + 2 \ |_i, 

2-^^r{£ + 2k + i)r{k)r{d + 2k)^ UJ-' k'^ + i)^ [ 2 J -' 1^ 2 j 
r{2k)r{l + k)r{£ + d + 2k + i) r (^ + fc) r (^±f±i) r (^±f±2) 

2^-^r(^ + 2fc + i)r(d + 2fc) r (|) 2^-^'=0Fr(2fc)2-^"''-^''y¥r(^ + d + 2fc + 1) 



2^-1^/rf 



r{2k)r{£ + d + 2fc + 1) 21-^^-2'= 0Fr(d + 2fc) 2-^-2^0^ j^(^ + 2fc + 1) 

The result for /g°° Gq(z) z'*"^ dz follows from Equation (3.326.1) in [9]. 
We are now ready to state the main result. 



Theorem 4.3 Let a be a positive real constant, k be a non-negative integer or half- 
integer and ^lii_k be the normalised equal area Wendland functions defined by (|25p with 
£ given by (|2.1|) . Then 

lim i^iAy) = Ga{y) (36) 

A;— >oo 



With the convergence being uniform for y € [0,oo). 
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Proof It follows from ([7]) that for a fixed y 

, d f°° / ^ \ d 

^ ^TdT / IJ'di'iA^) ~ Go.{z)\z''-' dz, (37) 

1 (2) Jo 

where we have used the following bound on the Bessel function [5l 10.14.4] 



Since the right hand side of ([37|l is independent of y, the result follows if we can 
show that 

/'OO 

Da,k-= I \^d^l,k{^)-Ga{z)\z''-^'^dz^O as fc^oo. 
Jo 

Now for arbitrary Z > we have, because J-d'4'i,k ^^nd Ga are non-negative, 

-Da,fc < / IJ'd^f.fcCz) -Gq(2)|z ^Mz+ / Td^i^k{z)z ~^dz+ / Ga(z)2 ^^c 



\Td^i,k{z)~Gc.[z)\z''-Uz+ / (Ga(^)-J-dV£,fcW)^''"'d2 



+2 / Ga(z)2''"Mz 

Jz 



IZ 
f-Z f-oo 

<2 \Tai^i^k{z)-Ga{z)\z'^-^dz + 2 Ga{z)z'^-^dz, (38) 

Jo Jz 

where in the second step we have used Lemma 14.21 

Given an arbitrary e > 0, we now choose Z sufHciently large to ensure that 

/ Ga{z) Z ^ dz 

Jz 



1.<1. 



For the first term in p8|l . we note that from Theorem l4.1l the integrand converges 
to zero uniformly for 2 in a bounded interval, thus there exists fci £ N such that 



z 



IJ'di^t^kiz) - Ga{z)\z'^-Uz < - W k>ki 
^ 



and hence 

for all k > ki, which completes the proof. 
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An interpretation of Theorem 14.31 in terms of probability distributions is that 

lim -T^hk f , f 1 . ) = -7^ exp (^^] , y > (39) 



k^ooy/2-Ka \^i.k[2t^) I V^Tva V ^'^ 

where the hmit on the right hand side is the right half of the Gaussian probability 
density with mean and variance a^. 
We note that since 

SiA^)-ixl^, (40) 

from (|34p . we could have used the right hand side of this expression instead of 
^i,k{<^) to define a change of variable to study the limiting case. 



5 Numerical results 



In this section we present numerical results regarding the differences between the 
appropriately scaled Wendland functions and the Gaussian limit established in 
Theorem l4.3l We also consider an interpolation example using both the normalised 
Wendland functions 4>£i; and the normalised equal area Wendland functions tpi,k- 



5.1 Difference with the limiting Gaussian 

Let the differences between the normalised equal area Wendland functions and the 
limiting Gaussian be 

E£,k{y) ■= i'e,k{y) ~ exp(-ay^) 



and let 



^i.k :=sup|£;f.fe(y)| 



Note that the change of variable used to define '(/'£, fc depends on the parameter a. 
Figure [3] shows plots of Eg )^{y) with a = 1. The upper plots are for d = 2 and 
k = 1.5 and k = 5.5 and the lower plots are for d = 3 and fc = 2 and fc = 6. 

In the absence of theoretical rates of convergence, we show numerical results. 
Figure [H shows e^ j. with q = 1 for k = 1, . . . , 50 and d = 3, 5, 7 and 9 for the 
original Wendland functions. Figure [5] shows egj. with a = 1 for k = 0.5, . . . ,49.5 
and d = 2, 4, 6 and 8 for the missing Wendland functions. Since a is just a scaling 
factor, the results do not vary in an essential way for different values of a. 

In all cases, we see rapid convergence of e^ j. to zero as the smoothness param- 
eter k increases. This is consistent with the theoretical convergence results. Note 
that e^ fc is not monotonically decreasing in k. We also remark that e^ j. is reached 
at different values of y as A; increases. 
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'^3,1 !<->■) 




£,,,.,(>') 




(a) d = 2,k = 1.5 



(h) d = 2,k = 5.5 





{c) d = 3,k 



(d) d = 3,fc: 



Fig. 3: Eg i^{y) with a = 1 and < j/ < 5^^a:(1). Subplots (a) and (b) are for the 
missing Wendland functions and subplots (c) and (d) are for the original Wendland 
functions. 



5.2 An interpolation example 



We consider an example, in which we show results obtained with both the nor- 
malised Wendland functions 0^ j, and the normalised equal area Wendland func- 
tions -(/)£ fe for different values of k. The aim of the example is to approximate the 
2-dimensional Franke function [6l p. 20] on [0,5]^. For k = 1,...,5 we consider 
interpolation as in dTJ and ([2]), using the normalised Wendland functions (f>ik s-nd 
the normalised equal area Wendland functions ■4'i,k with a = 2. We use 9x9 
and 17 x 17 equally spaced grids as the centres. The number of centres is thus 
N = 81 and N = 289. The L2 error was estimated using Gaussian quadrature with 
a 120 X 120 tensor product grid of Gauss-Legendre points and the Loo error was 
estimated by using a 360 x 360 equally spaced grid. Table |3] shows the L2 and Loo 
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0.005- X 



(a) d = 3 



(b) d = 5 



10 20 30 40 50 



(c) d = 7 



10 20 30 40 50 



(d) d = 9 



Fig. 4: e^ fc on a logarithmic scale with a = I, k = 1, . . . ,50 and d = 3, 5, 7, 9 for 
the original Wendland functions. 



errors, as well as the 2-norm condition numbers of the interpolation matrices. We 
also show the results with the limiting Gaussian of exp(— 2 y^), denoted by fc = oo. 

We see from the right-hand part of Table [3] that once the argument is properly 
scaled to give approximately constant effective support, increasing the smoothness 
has remarkably little effect on the error. On the other hand the condition number 
increases rapidly as the smoothness increases and is very large for the Gaussian 
limit. Taken together, these observations suggest that any benefit gained from the 
higher smoothness is likely to be offset by the increased condition numbers of the 
matrices. 



The results with the normalised Wendland functions (f>£k are in the left-hand 
part of Table[3l We can see that the condition number is decreasing as k increases, 
which is due to the decreasing magnitude of the non-zero elements away from the 
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9.5 19.5 29-5 39.5 49.5 

k 



(a) d = 2 



9.5 19.5 29.5 39.5 49.5 

k 



(b) d = 4 



9.5 19.5 29.5 39.5 49.5 



(c) d = 6 



9.5 19.5 29.5 39.5 49.5 



(d) d = 8 



Fig. 5: e^ fc on a logarithmic scale with a = 1, fc = 0.5, . . . ,49.5 and d = 2,4,6, 
for the missing Wendland functions. 



diagonal. This is due to the fact that as k increases the normalised Wendland 
functions <^f j, decay more rapidly with respect to r, as illustrated in Figure [T] 



6 Conclusion 



Compactly supported radial basis functions have proved very popular in scat- 
tered data approximation due to the resulting sparsity of the interpolation matrix. 
This paper has shown that for both the original and missing Wendland functions, 
the limit as the smoothness parameter goes to infinity, after suitable scaling and 
linear transformation, is a Gaussian RBF. 
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RBF: 4>e^k 


RBF: V^,fe 1 


N 


k 


L2 error 


Loo error 


K 


L2 error 


Lao error 


K 


'^min 


'^max 


81 


1 


2.25e-l 


6.96e-l 


1.71 


1.89C-1 


5.89e-l 


1.76cl 


1.55e-l 


2.74 




2 


2.6I0-I 


7.95e-l 


1.22 


1.86C-1 


5.78C-1 


3.14cl 


9.62C-2 


3.02 




3 


3.00e-l 


8.90C-1 


1.07 


1.87e-l 


5.79C-1 


4.96el 


6.50e-2 


3.22 




4 


3.36e-l 


9.73e-l 


1.02 


1.87e-l 


5.80C-1 


5.56cl 


5.98e-2 


3.30 




5 


3.63e-l 


1.03 


1.01 


1.87c- 1 


5.81C-1 


6.37cl 


5.29C-2 


3.37 




CO 






1.89C-1 


5.89C-1 


9.40cl 


4.03e-2 


3.78 


289 


1 


7.75e-2 


2.26e-l 


4.87el 


8.20e-2 


2.45e-l 


4.77e2 


2.25e-2 


10.74 




2 


7.48e-2 


2.00e-l 


5.64el 


7.98e-2 


2.11e-l 


3.89e3 


3.08e-3 


11.97 




3 


7.47e-2 


1.98e-l 


3.00el 


7.88e-2 


2.00e-l 


1.23e4 


1.04e-3 


12.78 




4 


7.59e-2 


2.12e-l 


1.24el 


7.71e-2 


2.04e-l 


9.88c4 


1.33e-4 


13.09 




5 


7.75e-2 


2.33e-l 


6.85 


7.61e-2 


2.09e-l 


3.77c5 


3.55e-5 


13.39 




oo 




7.23e-2 


1.90e-l 


2.13e9 


7.04e-9 


15.00 



Table 3: Results from the example in Section 15.21 showing L2 and Loo errors, 2- 
norm condition numbers k and minimum and maximum eigenvalues (Amin and 
Amax) when using the normalised Wendland RBFs 0^ j, and the normalised equal 
area Wendland RBFs 4'£.k with a = 2. 



In Figure [T] we saw that the (original) normalised Wendland functions exhibit 
faster decay with respect to r as the smoothness parameter k increases. This sug- 
gests the need for a change of variable, not only to have a well-defined limit as 
considered in this paper, but perhaps also in practical applications. Without a 
change of variable, in the case of interpolation we could have a nearly diagonal in- 
terpolation matrix and consequently high errors between the interpolation points. 

The number of interpolation points that fall within the support of an RBF 
is also related to the stationary approach to interpolation (e.g. [H] Section 12.2]), 
in which the goal is to keep the number of points in the support of each RBF 
approximately equal across different sets of centres. However here it is not the 
centres that are changing, but rather the basis functions that change with k. 

In Figure [2] we saw the (original) normalised equal area Wendland functions, 
whilst formally having support [0, <5f.fc(a)] that is different for each k, appear nearly 
identical. As a result, if the user wishes to compare the results using Wendland 
functions of different smoothness, the normalised equal area Wendland functions 
may be a more appropriate choice of RBF. 

We saw in Section 15.21 that whilst the normalised equal area Wendland func- 
tions with different k might appear comparable, and give similar accuracy, the 
increasing support as k increases causes decreased sparsity of the interpolation 
matrix, and consequently an increased condition number of the linear system. 
This leads us to conclude that there may be little benefit from considering high 
values of the smoothness parameter k. 



The results in the paper have illustrated the trade-off between approximation 
power and the condition number of the resulting linear system with Wendland 
functions of different smoothness. The issue of appropriate scaling and the selection 
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of a smoothness parameter when using the Wendland functions remains a complex 
issue in practice. 
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